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We investigate stability and dynamics of higher-order soli- 
tary waves in quadratic media, which have a central peak 
and one or more surrounding rings. We show existence of two 
qualitatively different behaviours. For positive phase mis- 
match the rings break up into filaments which move radially 
to initial ring. For sufficient negative mismatches rings are 
found to coalesce with central peak, forming a single oscillat- 
ing filament. 



Stability of optical solitary waves ('solitons') is one 
of the most important questions of theoretical nonlinear 
optics because of its direct connection with the possi- 
bility of experimental observation of solitons. Stability 
of the solitons in fully integrable systems naturally fol- 
lows from integrability. Solitons of the one-dimensional 
(ID) Nonlinear Schroedinger Equation (NLS), describing 
propagation of short pulses in a fibre with cubic nonlin- 
earity, are a well known example |1[. A wide range of 
the non-integrable hamiltonian models also have solitary 
solutions. For instance, equations describing parametric 
interaction in quadratic nonlinear media have solitary so- 
lution, which were pioneered in Refs. ^ and recently ex- 
plored in details from both theoretical and experimental 
sides because of their many interesting features, see e.g. 
p,0. In non-integrable systems the stability of ground- 
state solitary solutions is often governed by the deriva- 
tive of some integral invariant with respect to an associ- 
ated free parameter of the solution |Q-^. For example, it 
has been rigorously proved that for ground-state bright 
solitary solutions of the generalized NLS equation pos- 
itivity of the derivative of total energy with respect to 
the nonlinear wave number shift is a necessary and suffi- 
cient condition for stability B. Numerical and analytic 
studies indicate that this holds also for ground states in 
quadratic media [Q. 

Existence of higher-order solitary waves with bright 
and dark central spots surrounded by one or more 
rings was demonstrated for 2D NLS equation with pure 
Kerr Q and saturable gP-|lO| nonlinearities and also in 



quadratic nonlinear media | lOl- 13] . No universal stability 
criterion is known for higher-order bound states and their 
stability has to be studied individually in every case. It 
has been shown that for saturable nonlinearity higher- 
order bound states with bright and dark central spots 
are stable with respect to purely radial perturbations. 



obeying to criteria for ground states, but unstable with 
respect to azimuthally dependent perturbations, showing 
breakup of their rings into filaments |p[-p^ . Properties of 
solutions with dark central spot are strongly affected by 
their nonzero angular momentum and these properties 
are very similar for both saturable and quadratic nonlin- 
earities Q. 

Dynamics induced by instability of higher-order states 
is fascinating phenomena on its own and it is a natu- 
ral starting point for understanding pattern forming phe- 
nomena in the evolution of higher-order gaussian beams 
in nonlinear media |l3] . Primary interest of this Commu- 
nication is to address (for the first time to our knowledge) 
problem of stability and dynamics of higher order solu- 
tions with zero angular momentum in quadratic media. 
In particular we show that these solutions reveal new 
scenario of evolution which is absent for corresponding 
solutions in Kerr-like media. Namely for some parame- 
ters values symmetry-breaking instability leading to fil- 
amentation along rings replaces with novel symmetry- 
preserving instability resulting in coalescence of rings 
with a central peak. 

We consider interaction of first and second har- 
monic optical fields propagating in a dielectric medium 
with quadratic nonlinearity, under the conditions of 
type I phase matching and with negligible walk-off ef- 
fects. The corresponding hamiltonian is iJ = 
JJdxdy [i|V^i?i|2 + i|Vii?2|2 + p\E2\' - \{ElE'i + 

c.c.)], where Vj_ — idx + jdy and [3 is the normalized 
phase mismatch. All variables and parameters are dimen- 
sionless, and these scaled units are used throughout the 
text and in the figures. The evolution of the normalised 
field envelopes of the fundamental {Ei) and second (£'2) 
harmonics obeys the system of equations: 



idzEm = 



6H 

Je*: 



m = 1,2. 



(1) 



We look for non-diffracting solutions of Eqs. (1) in the 
form Em{z,x,y) = Am{r)e""-'^^ , where r = y^x^ + y^, 
Am are real functions and k is the nonlinear wave vec- 
tor shift. The existence condition of localised solutions 
with exponentially decaying tails is /t > max(0, —P/2). 
For any value of k in this range we were able to numeri- 
cally build higher-order many-ring solitary solutions with 
a bright central spot. Examples of spatial profiles of one- 
and two-ring solutions are presented in Fig. 1(a). For 
any finite number of rings, the fundamental field has ra- 
dial zeroes but the second harmonic field always remains 



positive, though having minima close to the zeros of the 
fundamental. In the limit /3 ^ 1 Eqs. (1) can be approx- 
imately reduced to an NLS equation for the fundamental 
field. Accordingly for increasing (3 the second harmonic 
field tends to carry less and less of the energy. The situ- 
ation is opposite for negative /3, when k values are close 
to the boundary of soliton existence. Dependences vs k 
of the energy invariant Q — J J dxdy (|i?ip -f 2|i?2p) are 
presented in Fig. 1(b). 

To study stability we consider 2D perturbations of 
these solutions in the general form Ejn{r,9,z) = e"""^ 
{A„,{r) + e+ (r)e^^+*''^ + e-*(r)e^*^-*''^). Here 9 is the 
polar angle and J is the azimuthal index of the pertur- 
bation. J must be an integer for azimuthal periodicity. 
Linearising Eqs. (1) and putting e^ = ?'''^'/j,„(?') gives 
the non-self-adjoint eigenvalue problem: 
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K, Lj2 — jRj - 2k - (3 and 
We seek unstable eigenmodes in 

the class of functions obeying ^^j^ = at r = and 
exponentially decaying at r — > cxd. Corresponding eigen- 
values belonging to the discrete spectrum can lie any- 
where in the complex plane outside the rays [i^c, ioo) 
and {^iQc, —ioo) which belong to the stable unbounded 
continuum, here flc = 7nin{K,2K + (3). Unstable modes 
have eigenvalues with ReXj > 0. They must always 
have a counterpart with ReXj < because of the hamil- 
tonian nature of our problem. Infinitesimal phase and 
translational transformations of the stationary solutions 
generate two eigenfunctions: /q^ — ±mAm for J = 
and /f^ = 1 dAm_ for J — 1. which are neutrally stable 
(Ao,i = 0). '' / 

We will mainly deal with solutions with only one ring 
outside the central peak. These show the main features 
of the dynamics of solutions with an arbitrary number 
of rings. First we consider symmetry-preserving pertur- 
bations, J = 0. Using asymptotic techniques developed 
for the ground-state solution B , it can be shown that the 
neutrally stable mode branches at the point d^Q = giv- 
ing instability for d^Q < 0, see Fig. 1(b). Thus we can 
conclude that the standard stability criterion for ground 
states y,y is also a necessary condition for stability of 
higher-order bound-states. This instability is related to 
the existence for c^^Q > of a pair of eigenmodes with 
purely imaginary eigenvalues (with opposite signs) lying 
in the gap (— i51c,i^c)- At the point d^Q = these 
eigenmodes coincide with the neutral mode and for more 
negative /3 appear again but with purely real eigenvalues 
of opposite sign, signifying instability. For the ground- 
state this is the only instability scenario and these dis- 
crete eigenmodes disappear into the continuum for large 
/3>0 0. 



In the present system, we have undertaken numerical 
investigation of eigenvalue problem (2). The case J — 
reveals two pairs of discrete eigenmodes. Interplay be- 
tween them leads to a new bifurcation scenario, which 
we study for different values of f3 for fixed k — 3. Chang- 
ing K at fixed (3 has no qualitative effect due to scaling 
properties of (1). However introducing this scaling mod- 
ifies the stability criterion di^Q > y, which we prefer 
to avoid. 

Real and imaginary parts of key eigenvalues from the 
discrete spectrum are plotted vs /? in Fig. 2. In the 
limit of large f3 we found one internal eigenmode (line 

1 in Fig. 2) but at /3 ~ 4.75 another internal eigen- 
mode (line 2) emerges from the continuum. On emer- 
gence mode 2 has ImXo = k = 3, but as /3 is decreased 
the eigenvalues of the two modes bring together, as Fig. 

2 shows. They fuse at /3 ~ —0.82 to form two pairs of 
eigenfunctions with complex conjugate eigenvalues, giv- 
ing onset of an instability (lines 3,4). At /? ~ —2.06 a 
reverse bifurcation takes place (lines 5,6). One eigen- 
mode (Fig. 3(a)) then has a purely imaginary eigenvalue 
(line 6) until it loses its stability at /3 ~ —5.61 (line 7) 
where d^Q = 0. This is the standard instability scenario 
described above. The other eigenmode undergoes a sim- 
ilar bifurcation, but at /3 ~ —2.17 (lines 5,8) which is 
well before the point d^Q — 0. At this bifurcation the 
eigenfunction profiles. Fig. 3(b), are quite different from 
those of the neutrally stable eigenmode. Thus, unlike the 
previously- known case (lines 6,7 and Fig. 3a) this new 
instability cannot be captured by asymptotic expansion 
around that neutral eigenmode. In both cases the unsta- 
ble eigenvalues reach a maximum then go steeply to zero 
near the existence limit of solitary solutions, (3 — —2k. 

The cascade of symmetry preserving bifurcations pre- 
sented in Fig. 2 is somewhat similar to that for TEi 
mode instability in a planar waveguide (ID geometry) 
with Kerr nonlinearity |15| where joint action of the re- 
fraction index discontinuties and field nodes leads to in- 
stability. In our situation the new symmetry-preserving 
instability develops in the region where the nodeless sec- 
ond harmonic starts to dominate over the fundamental, 
which has one or more nodes. 

For symmetry-breaking perturbations, J 7^ 0, the sta- 
bility properties are not so rich as for J = 0. For J — 1 
our numerics reveals the presence of a neutral mode and 
a pair of eigenmodes in the discrete spectrum with purely 
real eigenvalues, one of which is responsible for instabil- 
ity. We did not find any exchange of stability of these 
modes. For every J from 2 to 5 we find such a pair of 
discrete eigenmodes with purely real eigenvalues and all 
modes for J > 5 belong to the continuum. The unstable 
perturbations for J = 3, 4 are localised around the ring 
of the bound-state in a manner similar to what happens 
in saturable media |g]. 

To show how the character of instability of the one- 
ring solution depends on phase mismatch parameter we 
plot in Fig. 4 growth rates vs f3 for all unstable eigen- 
modes. For phase mismatches from the cascading limit 



down to /3 '^ —3 symmetry breaking instabilities with 
J = 3, 4 are dominant. However, sufficiently far from 
the NLS limit our new scenario, with azimuthally homo- 
geneous perturbations dominant, is realised. We stress 
again (see discussion above) that this symmetry preserv- 
ing instability is not related to violation of the criterion 
dttQ > 0. Note that in the limit /3 ;2> 1 the J — internal 
eigenmode exists and the J = 1 eigenmode has non-zero 
growth rate (Fig. 4(b)). 

One expects the propagation dynamics of solitary 
states to be mainly determined by the perturbation 
eigenmode with maximal growth rate. To examine this 
we performed an extensive series of numerical simula- 
tions of Eqs. (1), using both polar and cartesian grids. 
Predictions based on our stability analysis are in good 
agreement with the results of our simulations. An exam- 
ple of noise-stimulated break-up of a one-ring solution 
into three filaments is shown in Fig. 5(a). We plot the 
real part of the fundamental field profile rather than the 
intensity distributions to show that the daughter solitons 
formed from the ring are out of phase with the central 
one. Radiation losses in the break-up are quite small, 
so that the initial energy Q is mostly divided among the 
daughter solitons. Their diameters are comparable to the 
width of the initial ring. For /3 values where growth rates 
for J — 3 and J = 4 are almost equal the simulation re- 
sults depended on the particular noise realisation, but we 
mostly observed the ring forming four filaments, one of 
which was usually less intense than the others. 

Throughout the whole range of parameters where the 
symmetry-breaking instability is dominant we observed 
repulsion between the central spot and daughter fila- 
ments, which results from the fact they are out of phase 
|l6| , see Fig. 5(a). This repulsive force makes the outer 
filaments move out along radii (Fig. 5(b)), in contrast to 
the tangential motion of daughter solitons after breakup 
of one-ring solitary waves carrying non-zero orbital an- 
gular momentum llG] , where inter-soliton forces are neg- 
ligible in comparison to the need to conserve angular mo- 
mentum. 

Our stability analysis predicts a novel symmetry- 
preserving instability scenario where the J = eigen- 
mode dominates. This prediction is indeed confirmed by 
the simulations. For example, at /3 = —4.2, instead of 
fragmentation we observed coalescence of the ring with 
the central spot to form a single filament. After transient 
dynamics this filament forms an oscillating solitary wave, 
see Fig. 6. These undamped pulsations are related to the 
existence of an internal eigenmode of the ground-state 
solution Jm). 

Considering now two-ring solitary solutions, we present 
the growth rates for the dominant eigenmodes and an ex- 
ample of symmetry- breaking instability, see Fig. 7. Gen- 
eral features of the dynamics are qualitatively similar to 
the one-ring situation. 

The evolution of filaments following a symmetry break- 
ing instability of peak-and-ring solitary solutions in sat- 
urable Kerr media |q] is a question which has not previ- 



ously been examined. In simulations of this problem we 
observed the same sort of dynamics as described above 
for quadratic media, but with no coalescence phenomena. 

In summary, we have undertaken a detailed analysis of 
stability of cylindrically-symmetric higher-order solitary 
waves due to parametric interaction in quadratic non- 
linear media. For a wide range of positive mismatches, 
symmetry-breaking instability of the rings is predicted, 
and confirmed by simulations which show that the in- 
stability leads to filamentation into daughter solitons 
which are repelled radially from the central spot. For 
sufficiently negative phase mismatches we predict that a 
new symmetry- preserving instability becomes dominant. 
This is confirmed in simulations, in which the rings are 
found to coalesce with the central filament, forming an os- 
cillating, single-peaked, solitary wave. To our knowledge 
this is the first explicit example of a symmetry-preserving 
instability of 2D self-trapped beams in bulk media differ- 
ent from the Vakhitov-Kolokolov scenario |5| . 
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FIG. 1. (a) Radial profiles of one- and two-ring solitary 
waves. Full (dashed) lines are for Ai (^2). (b) Total energy vs 
K for one-ring (full lines) and two- ring (dashed lines) solitary 
waves. The negative values of /3 are chosen so as to give 
OkQ = at k = 3. 




FIG. 2. Real (full lines) and imaginary (dotted lines) parts 
of the eigenvalues of J = eigenmodes vs P, k — 3. 



U.2U 


% 


0.3 





0.15 


\ (a) 


0.2 


V (b) 


0.10 




0.1 


i 


0.05 


^\ V 


0.0 


%2\., 


0.00 
0.05 




i^;,.-^ 


-0.1 







50 100 150 200 250 300 



FIG. 4. Growth rates of the maximally unstable eigen- 
modes vs /3 for one-ring solitary solutions. Dotted line in 
(b) displays ImXo of the J = internal eigenmode marked 
by Line 1 in Fig. 2. k — 3. 



14 




12 


(a) - 


10 


• 


8 


• 


6 


• 


4 




2 

n 





2 4 6 8 101214 




10 15 20 25 



FIG. 5. (a) Real part of the fundamental harmonic field at 
a late stage of a simulation of the symmetry-breaking process, 
(b) Superimposed images of its transverse intensity distribu- 
tion at different values of z, showing radial trajectories of the 
daughter solitons. I3 — —1, k — 3. Brightness and size of 
central spot in Fig. (b) are exaggerated by the superposition 
of multiple images. 
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FIG. 6. (a) Initial (dashed lines) and post-coalescence 
(solid lines) radial profiles of fundamental and second har- 
monics simulated where the symmetry-preserving instability 
is predicted to dominate, (b) Corresponding evolution of 
\Ei^2\ at r = vs propagation coordinate z. n = 2>, (3 = —4.2 
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FIG. 3. (a) The internal eigenfunction corresponding to 
branch 6 in Fig. 2 at the point /3 = —5.58, slightly before the 
bifurcation point d^Q = 0. (b) The internal eigenfunction 
corresponding to the bifurcation point at /3 = —2.17 where 
branches 5,8 in Fig. 2 meet. Dots mark the neutrally stable 
eigenmode. k = 3. 
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FIG. 7. (a) Growth rates of the dominant unstable eigen- 
modes vs /3 for two-ring solitary solutions, (b) Real part of 
the fundamental field at a late stage of a simulation of the 
symmetry-breaking process. l3 — 1, k — 3. 



